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Abstract 

The London ground-state energy formula as a function of number density p for a system 
of boson hard spheres of diameter c at zero temperature (corrected for the reduced mass 
of a pair of particles in a "sphere-of-influence" picture) generalized to describe fermion 
hard-sphere systems with four and two intrinsic degrees of freedom such as 3 He or neu- 
tron matter and symmetric nuclear matter, respectively, is proposed as the crystalline 
energy branch for hard-sphere systems. For the fluid branch we use the well-known, ex- 
act, low-density equation-of-state expansions for many-boson and many-fermion systems, 
appropriately extrapolated to physical densities. Here, via a double-tangent construction 
the crystallization and melting densities for boson and fermion hard spheres are deter- 
mined. They agree well with variational Monte Carlo, density-functional, and Green 
Function Monte Carlo calculations. 

PACS: 05.30.-d; 21.65.+f; 67.90.+Z 

Key words: boson and fermion hard-sphere systems; nuclear and neutron matter; freezing 
transition. 
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1 Introduction 



As in the classical hard- sphere system where a "freezing" transition from a fluid to a crystalline 
phase was found [1] to occur as one increased density, one expects the same phenomenon to 
occur in quantum systems, whether of bosons or fermions. 

An analytical formula for the ground state energy E of a boson hard-sphere system for all 
densities p was proposed by London [2] as 

^ /Ar 2nh 2 c 11 

EN= -p: TT ~, 1 

m ( p -i/3_ /9o -V3 )2(p - 1/3 + fep -i/3 ) 

where m is the particle mass, c is the hard-sphere diameter, with p = a/2/c 3 the ultimate 
density [3] at which a system of identical classical hard spheres close packs in some primitive- 
hexagonal arrangement, e.g., face-centered-cubic, and b = (2 5 ^ 2 /n) — 1. The justification given 
for Eq. (1) is that it reduces smoothly to well-known limiting expressions at both low and high 
densities, namely 

E/N ► (2nh 2 /m)pc, (2) 

E/N — + A(h 2 /2m)(p- l ^-p 1/3 )- 2 . (3) 

Here A = 7r 2 /2 1//3 ~ 7.834 is a constant called the residue of the pole at close packing. The 
low-density leading term (2) is the celebrated Lenz [4] term, calculated by him as the leading 
correction to the energy due to an "excluded volume" effect. The limiting case Eq. (3) comes 
from the lowest Schrodinger equation eigenvalue of a particle in the spherical cavity. It is just 
the kinetic energy of a point particle of mass m inside a spherical cavity of radius r — c, where r 
is the average separation between two neighboring hard spheres and r = (V2/ p) 1 ^ 3 by assuming 
a primitive-hexagonal arrangement for the cavities. 

It was found [5] that the derivation of the high-density extreme of the original [2] (boson) 
London equation (1) contains one fundamental error, the spherical cavity of radius r — c alluded 
to above in reality refers to a "sphere of influence" of two particles. Thus, the particle mass 
used in obtaining (3) should refer to the reduced mass m/2. This yields the constant 

b = 2 3 / 2 /tt - 1 (4) 

instead of 2 5 / 2 /tt - 1 given by London (1). The result (1) with (4) is designated the modified 
London (ML) equation, which continues to satisfy (2) as this is independent of the constant b. 
The residue A in (3) then becomes 2 2 / 3 7r 2 ~ 15.7 and fully agrees with the empirical residue of 
15.7 ± 0.6, extracted by Cole [6] from high-pressure crystalline-branch data in 3 He, 4 He, H2 and 
D 2 systems. Moreover, this ML equation agrees dramatically better than the original London 
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(L) equation with Green Function Monte Carlo (GFMC) [7] computer-simulation datapoints 
for both fluid and crystalline branches of the boson hard-sphere system. 

A generalized London equation has also been proposed [8] for N-fermion hard-sphere sys- 
tems, with v intrinsic degrees of freedom for each fermion. Two differences appear with respect 
to the boson London formula: a) unlike the boson case, the ground state kinetic energy is 
nonzero and is added as a well-known [9] //-dependent leading term, and b) the constant b is al- 
lowed to be //-dependent, namely b(u) = [{v — l)/v](b + 1) — l.This form ensure a v- independent 
energy at close-packing where, since the spheres can be labelled indistinguishability and thus 
particle statistics are lost as expected in this classical limit. Replacing the value of the constant 
b by the new one (4) gives a generalized modified London equation (ML„) 

n/ „ r _ 2/3 (V-l\2lxh 2 C 1 1 

E N = C u p 2/3 + [ m TTT- 5 

with 

C^^f-^o (6) 
10m V v J 

where in the limit v — > oo the new constant b(u) — > b and (5), i.e., goes over into the boson 
case stated above. The low-density limit of Eq. (5) becomes 

E/N-fC.?» + (!^)™?.pc (7) 
p-^o V v J m 

where the second term on the rhs is the Lenz term for //-component fermions. On the other 
hand, for p — > p = \[2/c 3 one sees that (5) reduces to (3) as it should. In other words, 
hard-sphere fermions, bosons or "boltzons" must close pack at the same density. From this it 
follows that the residue for bosons or fermions is the same and equal to 2 2 ^ 3 tc 2 ~ 15.7, in line 
with Ref. [6]. 

In addition to the Lenz term (2) for the low-density fluid branch, several higher-order 
corrections to the ground state energy per particle have been derived using quantum field- 
theoretic many-boson perturbation theory [10], namely 

E/N = ^V(i + Cl {pc 3 ) 1 ' 2 + C 2 pc 3 ln(pc 3 ) + C 3 pc 3 + o(pc 3 )), (8) 

where d = 128/15^, C 2 = 8(4^/3-^), but C 3 is an as yet unknown constant. Here c is the 
S-wave scattering length in general, and becomes the sphere diameter for a hard-core potential. 
The series is clearly not a power series expansion, and is at best an asymptotic series. 

Similarly, for an iV-fermion hard-sphere system the corresponding series is [11] 

E/N = + dikpc) + C 2 (k F c) 2 

+ [C 3 ro/2c + C 4 Ai(0)/c 3 + C 5 ](k F c) 3 + C 6 (k F c) 4 \n(k F c) 

+ [C 7 r /2c + C 8 A'^0)/c 3 + C 9 ](k F c) 4 + o(k F c)y (9) 
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where the Cj (j = 1, 2, 9) are dimensionless coefficients depending on v and are given in Ref. 
[12] for v = 2 and v — 4. The Fermi momentum hkp is defined through the particle number 
density p = N/Q — vkp/6ir 2 , with f2 the system volume. 

Unfortunately, both low-density expansions (8) and (9) lack accuracy at moderate to high 
densities, including the saturation (or equilibrium, zero-pressure) densities of liquid He [y = 
oo) and liquid 3 He (is = 2) or nuclear matter {y — 4). However, one can extrapolate the series 
for hard-sphere systems to physical and even to close-packing densities through the use of Pade 
and so-called "tailing" [13] approximants. This method, called Quantum Thermodynamic (or 
Van der Waals) Perturbation Theory (QTPT) [14], has provided fairly accurate representations 
of the fluid branch of the equation of state [15] beyond the presumable phase transition densities 
but without enough credibility as one approaches close packing. This is clear since one does not 
possess a ground-state energy function with implicit information of both fluid and crystalline 
branches. Thus, for the crystalline branch we employ the generalized modified London equations 
(1) which reproduce the correct closest packing density value p = a/2/c 3 as well as good behavior 
near closest packing as suggested by a correct residue value. For the fluid branch we use the 
Pade extrapolations based on the series (8) and (9). 

This paper unfolds as follows. In Sections 2 and 3 we construct the fluid branches for hard- 
sphere bosons and fermions, respectively. In Sec. 4 we use them together with their respective 
modified London equation crystalline branches to determine, via a double-tangent construction 
ensuring equality of pressure in both phases, the melting and crystallization densities as well 
as energy and density changes. In Sec. 5 we state our conclusions. 

2 Boson hard-sphere fluid branch 

In order to extrapolate to higher densities the low-density series (8), we define 

E /N = —pce (x), x=(pc 3 ) 1/2 , (10) 
m 

with 

e (x) = 1 + dx + C 2 x 2 Inx 2 + C 3 x 2 + 0(x 3 Inx 2 ). (11) 
Alternatively, one can analyze the series 

e 1/2 (x) = 1 + F 1 x + F 2 x 2 Inx 2 + F 3 x 2 + 0{x 3 Inx 2 ) (12) 

for x <C 1, where the F's are expressible in terms of the C's, with C 3 and consequently F 3 
unknown. Values of the C's and F's are given in Table 1. We analyze eQ 1 ^ 2 (x) instead of 
eo(x) to ensure that zeros in its extrapolants e (x) are second-order poles in energy as we 
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expect at close packing (CP). Then we examine the twelve extrapolants [16] from the series 
(12) with three terms. Fitting extrapolants to go through the four GFMC data points [7] 
allows estimating a good value for the last coefficient F 3 . Forms VII and XI both satisfied these 
conditions but the mean square deviation with respect to the four GFMC data points was least 
with form XI. Thus we use it as our best extrapolant. So the ground state energy for boson 
hard spheres is represented (symbol =) by 

E/N=—pce (x), (13) 
m 

with 

-i/2 YT , v 1 + F 2 x 2 \nx 2 

60 =XI{X)= 1-F 1X -(F 3 -F?) X 1 (14) 
and F 3 = —27.956. In Fig. 1 this is plotted as a full curve and labeled XI for the fluid branch. 
Open circles and squares are GFMC data for fluid and crystalline branches, respectively. The 
dashed curve is the Modified London formula (1) representing the crystalline branch. Dots are 
Diffusion Monte Carlo (DMC) [17] calculations spanning a wider range of densities in the fluid 
region than the GFMC data. From DMC data we see that although our expression for the 
fluid branch could be better for intermediate densities, it agrees well with DMC and GFMC 

— 1/2 

calculations around the freezing transition. This e predicts a random close packing (RCP) 
density p/p = 0.776 which is only ten percent below classical RCP p/p ^ 0.86 [3] expected to 
be the ultimate CP density also for quantum hard sphere fluids just from the fact that particles 
at CP are localized so that indistinguishability does not hold and thus statistics are irrelevant. 



3 Fermion hard-sphere fluid branch 

For fermion hard-sphere systems we write (9) as 

3 h 2 k 2 

E/N = -— -^eo(x); x = k F c, p = N/tt = isk 3 F /6ir 2 , (15) 



Table 1. For bosons, Cj and Fi coefficients appearing in Eqs. (11) and (12), 
respectively. Numbers in quotation marks are adjusted as indicated in text. 



Bosons [y = oo) 


i =1 


2 


3 


Cj 


4.81441778 


19.65391518 


"73.296" 


Fi 


-2.40720889 


-9.826957589 


"-27.956" 
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0.0 0.1 0.2 0.3 0.4 

P/Po 

Figure 1: The quantity £o 1/2 = ^27ch 2 pcN/mE = [1 - (p/p ) 1/3 ]^l + b (p/p ) 1/3 function 
of p/po for boson hard sphere systems: XI is the fluid branch (14), ML is the Modified London 
formula (1), open circles and squares are GFMC data for the fluid and crystalline branches 
respectively, and dots are DMC calculations. 

with 

e (x) = 1 + dx + C 2 x 2 + (C 3 /3 + C4/3 + C 5 ) x 3 

+ C 6 x 4 lnx+{C 7 /3-C 8 /3 + C 9 )x 4 + o{x 4 ). (16) 

For v = 2 . C 6 = [11] so that (16) becomes 

e (x) ~ 1 + D 1 x + D 2 x 2 + D 3 x 3 + D A x 4 + o(x 4 ) (17) 

for x < 1. With D's in terms of C's and p = kp/37i 2 . As in the boson case, instead of eo we 
write 

e 1/2 ~ 1 + Fi x + F 2 x 2 + F 3 x 3 + F A x 4 + F 5 x 5 + o(x 5 ), (18) 

where the Fj's depend algebraically on the D^s in a simple manner, F 5 being unknown. Values 
of Di and Fj are given in Table 2. This series is a simple power series and we use it to construct 
the usual Pade extrapolants. The approximants to (18) with four terms beside the trivial 
one were analyzed in Ref. [18] concluding that the best was the Pade [4/0] (x). However, 
this function does not have a zero in the region of physical interest, i.e. < p/po < 1 or 
< x < 3.47, and so the energy does not manifest a CP as it should. It was necessary to 
introduce the fifth term in (18). Its five two-point Pade approximants [3//2](x) were analyzed 
and F 5 adjusted to ensure a zero. The position of the zero and the approximant were chosen 
in such way that the QTPT applied to 3 He with the Aziz inter-atomic potential reproduces 
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the corresponding GFMC [20] data. Eventually, the best extrapolant was the two-point Pade 
approximant [19] 

6o 1/2 = [3//2](x). (19) 

Hence E/N becomes 

E/N = ^([3//2](x))- 2 , (20) 

with a CP density pj p = 0.732. The coefficient F 5 is listed Table 2 between quotation marks. In 
Fig. 2 e 1/2 = (3/i 2 (67r 2 p/^) 2/3 ^/10m J B) 1 /2 = i + ( 2 0tt(z/ - l)/3z/)(2 1 / 4 z//67r 2 ) 2 / 3 {[(p/p )- 1 /3 

l] 2 [(p/p ) -1 ^ 3 — K^Kp/Po) 2 / 3 } -1 as a function of p/p is plotted for fermion hard sphere. For 
v = 2 the fluid branch [3//2] (full curve) is close to the numerically exact Ladder (Ladder) [21] 
(open squares), Variational Fermion-Hypernetted-Chain (VFHNC) [22] (plus sign marks) and 
L-expansion data [23] (open triangles). Fig. 3 is a enlargement of Fig. 2 where we show this 
excellent agreement over the whole range of available data. 




P/Po 

Figure 2: The quantity £q 1/2 = {^{^p/vf^N/lQmE) 1 / 2 = 1 + 20vr(z/ - l)/3u x 
(2 1 /V/67T 2 ) 2 / 3 {[(p/p )" 1/3 - l] 2 [(p/po) _1/3 -fe(^)](p/po) 2/3 } _1 as a function of p/p for fermion 
hard sphere with v — 4 (XII) and v = 2 ([3//2]), are full lines. Dashed lines are the corre- 
sponding Modified London formulas. 

For v = 4 . Eq. (16) becomes 

e ( x ) = 1 + D lX + D 2 x 2 + D 3 x 3 + D 4 x 4 Inx + D 5 x A + o(x 5 ) (21) 
for x 1. By analogy with bosons and fermions {v — 2), we analyze 

e 1/2 = 1 + F l x + F 2 x 2 + F 3 x 3 + F A x 4 lnx + F 5 x 4 + o(x 5 ) (22) 



7 



Table 2. For v — 2, Dj and Fj coefficients appearing in Eqs. (17) and (18), 
respectively. Numbers in quotation marks are adjusted as is indicated in text. 



v = 2 


i =1 


2 


3 


4 


5 


A 


0.353678 


0.185537 


0.384145 


-0.024700 


"-0.265544" 


F% 


-0.176833 


-0.045863 


-0.156677 


0.109672 


"0.130830" 



with all stated Fj known. Values of Di and Fj are given in Table 3. Unlike the v = 2 case, 
this series is not a pure power series as it contains logarithmic terms. Its tailing approximants 
are giving in Table III of Ref . [13] . We choose form XII to avoid approximants with spurious 
unphysical poles within the interval of physical densities, as well as to avoid residues falling 
outside the rigorous interval 

1.63 < A < 27. (23) 

This interval was obtained [24] for regular CP (face-centered-cubic or hexagonal-close-packing) 
by generalizing the exact calculation for a simple cubic lattice based on three mutually perpen- 
dicular linear lattice which gives A = ir 2 . 



1.0- 

CM 
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E 

O 0.8- 



cm 



CM 

CD 



0.6 



CM 
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CO 0.4 



v = 4 \\v 
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[3//2] 
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\ ML 2 


+ VFHNC(v = 2, 4) 
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a L-expansion (v = 2, 4) 

I I I I I 
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\\ \ 

l 
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Figure 3: Enlargement of Fig. 2 at low densities. 
Hence, E/N can be written 



with 



-V2, 



x) = Xll(x) 



1 + (Fi - F 3 /F 2 )x + (F 2 - F 1 F 3 /F 2 )x 2 
1- (F 3 /F 2 )x - F iX 4 \nx 



(24) 
(25) 
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Table 3. For v — 4, Dj and F { coefficients appearing in Eqs. (21) and (22), 
respectively. 



v = 4 


i =1 


2 


3 


4 


A 


1.061033 


0.556610 


1.300620 


-1.408598 


Fi 


-0.530517 


0.143867 


-0.5806558 


-0.704299 



This is plotted in Fig. 2 and 3 as a full curve and labeled XII. We also plot the corresponding 
VFHNC data (plus signs) and L-expansion data (open triangles). In terms of energy, ours are 
slightly below the VFHNC ones. The agreement gets better for lower densities. However, the 
XII approximant is just above the L-expansion data over the range of densities where they are 
available. 

4 Kirkwood phase transition 

Fluid-to-crystal (or freezing) phase transitions for boson hard-sphere systems have been de- 
termined with VMC, GFMC or Density Functional Theory (DFT) [25]. However, for fermion 
hard-sphere systems such calculations have not been done. The purpose of this paper is to 
provide very simple but hopefully accurate analytical expressions for the ground-state energy 
of quantum hard-sphere systems in general so as to reproduce crystallization transitions or to 
give the values for the parameter involved, as well as other thermodynamic properties such as 
pressure and compressibility. 

Here, via a double-tangent construction we locate the crystallization pi and melting p s 
number densities for bosons and for fermions (with v = 2 and 4). This process is schematically 
sketched in Fig. 4 for fermions with v = 2, where both fluid and crystalline energies per particle, 
in h 2 /mc 2 units, as a function of p/po, are shown close to the phase transition. Dots denote 
crystalline and fluid phase Maxwell-construction endpoints on the pressure-volume diagram, 
obtained by constructing a common tangent to both energy branches ensuring equality of 
pressure in both phases. 

For boson hard-spheres, in Table 4 we list our calculation (ML) together with the VMC, 
GFMC and DFT ones. All of them agree reasonably well with each other. Here A(pc 3 ) and 
A(E/N) signify the change of density and energy, respectively, between both phases. For 
fermion hard spheres with v — 4, e.g., symmetric nuclear matter, we find a crystallization 
density given by pi c 3 = 0.056 and a melting density given by p s c 3 = 0.069. This crystallization 



9 



Table 4. Crystallization and melting densities for boson hard sphere system. NA 
means not available. 



Bosons 


Pi c 3 


ps c 3 


A( Ps c 3 ) 


A(E/N) 


VMC [26] 


0.23 ± 0.02 


0.25 ± 0.02 


0.02 


1.28 


DFT [25] 


0.246 


0.284 


0.38 


2.668 


GFMC [7] 


0.25 ± 0.01 


0.27 ± 0.01 


0.02 


NA 


ML 


0.21 ± 0.02 


0.23 ± 0.02 


0.02 


1.37 



Table 5. Crystallization and melting densities for fermion [y = 4, 2) hard sphere 
system. 



Fermions 


Pic 3 


p s c 3 


A(pc 3 ) 


A(E/N) 


i/ = 4 


0.056 ± 0.005 


0.069 ± 0.006 


0.013 


0.29 


v = 2 


0.408 ± 0.001 


0.422 ± 0.001 


0.014 


1.589 



density is of the order of that of matter nuclear (0.32 ± 0.02 fm~ 3 ) [27, 28, 29] when we use 
the hard-sphere radius c = 0.64 fm determined through the equation of state of Jaqaman et 
al. [30] for nuclear matter interacting through a Skyrme effective interaction and adjusted to 
reproduce energy-per-nucleon and saturation density, i.e., -16 Mev and 0.17 fm -3 , respectively. 

For fermion hard spheres with v — 2, e.g., helium three or neutron matter, we find p\ c 3 = 
0.408 at crystallization and p s c 3 = 0.422 at melting. These density values agree well with the 
experimental (pi c 3 = 0.39, p s c 3 = 0.41) and calculated [31, 32] values for helium three. On 
the other hand, taking for the hard-core diameter a value of c = 0.4 fm, and assuming that 
an additional attractive potential would not change the crystallization density of real neutron 
matter too much, one could predict a crystallization density for neutron matter of about 6.38 
fm -3 (or mass density mp = 10.85 x 10 15 g-cm~ 3 , with m the neutron mass) which not only 
agrees with a prediction [33] that pure neutron matter crystallizes when its mass density exceeds 
1.5 x 10 15 g-cm~ 3 but also justifies Pandharipande's claim [34] p. 695, of no observation of a 
crystal phase below 4.2 fm~ 3 . In Table 5 we summarize our freezing transition parameters for 
fermion hard spheres. 
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0.26 0.28 0.30 0.32 

P/Po 

Figure 4: Energy E/N (in units of ft 2 /mc 2 ) as a function of p/po for fermion (y = 2) hard sphere 
system and Maxwell double tangent construction. Dots mean crystalline and fluid edges. 

5 Conclusions 

In summary, we have proposed analytical expressions for crystalline and fluid branches of both 
boson and fermion quantum hard-sphere systems. We have been able to calculate the fluid- 
crystal transition parameters for bosons and fermions such as symmetrical nuclear matter and 
neutron matter. The results support the controversial [34] possibility of crystallization in nu- 
clear matter as well as in neutron matter. 
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